module GaussJordan_module
	implicit none
contains
	subroutine	gj_gaussJordan(matrix, solution)
		real, allocatable, dimension(:,:), intent(in)	::	matrix
		real, allocatable, dimension(:), intent(out)	::	solution
		real, allocatable, dimension(:,:)				:: 	gauss_mat
		integer		::	rows, cols
		integer		::	i, j, index_of_pivot
		integer		::	m,n

		rows = size(matrix, dim=1)
		cols = size(matrix, dim=2)
		allocate(gauss_mat(rows, cols))
		gauss_mat = matrix
		!generating the diagonal form selecting the pivot
		do j = 1, cols-1
			index_of_pivot = absMaxOf(gauss_mat(j:rows,j))	!pivot is the absolute maximum of the jth column			
			index_of_pivot = index_of_pivot + j - 1
	
			call swapRows(gauss_mat, j, index_of_pivot)
			
			if( gauss_mat(j,j) == 0	) then
				allocate(solution(rows))
				solution = (/(0, i=1,rows)/)
				print*,"No solution"
				return
			end if		
			gauss_mat(j,:) = gauss_mat(j,:)/gauss_mat(j,j)
															
			do i = 1, rows
				if(i == j) cycle
				gauss_mat(i,:) = gauss_mat(i,:) - gauss_mat(i,j)*gauss_mat(j,:)
			end do
		end do
		allocate(solution(rows))
		solution = gauss_mat(:,cols)
	end subroutine gj_gaussJordan

	subroutine gj_checkError(matrix, solution)
		real, allocatable, dimension(:,:), intent(in)	:: matrix
		real, allocatable, dimension(:), intent(in)		:: solution
		real, allocatable, dimension(:)					:: error
		integer	:: i, j
		integer	:: rows, cols
		real	:: sum, mean

		rows = size(matrix, dim = 1)
		cols = size(matrix, dim = 2)
		allocate(error(rows))
		sum = 0
		do i = 1,rows
			error(i) = abs(DOT_PRODUCT(matrix(i,:cols-1),solution) - matrix(i,cols))
			sum = sum + error(i)
		end do
		mean = sum/real(rows)
		print*,'mean=',mean
	end subroutine gj_checkError

	subroutine	swapRows(matrix, m, n)
		real, allocatable, dimension(:,:), intent(inout)	::	matrix
		integer, intent(in)								:: 	m, n
		real, allocatable, dimension(:)		::	temp

		allocate(temp(size(matrix, dim=2)))
		temp 		= matrix(m,:)
		matrix(m,:) = matrix(n,:)
		matrix(n,:)	= temp
	end subroutine swapRows

	function absMaxOf(column)
		real, dimension(:), intent(in)	:: 	column
		real	::	temp(size(column))
		integer	::	absMaxOf
		integer	::	i
		real	::	max
		temp = abs(column)
		absMaxOf = maxloc(temp, 1)
	end function absMaxOf
end module  GaussJordan_module
